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We report progress for the development of MH4D for the first and second quarters of 
FY2004, December 29, 2002 - June 6, 2003. The present version of MH4D can now 
solve the full viscous and resistive MHD equations using either an explicit or a semi- 
implicit time advancement algorithm. This fulfills with what we promised in the 
proposal. 

In this report we describe progress in the following areas 

Presentation of a Poster at the EGS-AGU-EUG Joint Assembly in Nice, 
France, April 6-11, 2003. 

Presentation of a Poster at the 2003 International Sherwood Theory 
Conference in Corpus Christi, Texas, April 28-30 2003. 

Code Development. 

o Implementation of the MHD equations, 
o Implementation of the semi-implicit algorithm. 

Validation of the New Features Implemented in the Code 
o The Sod’s problem, 
o AlfV&i waves in a cylinder. 

Details are given in the following sections. 


I. PRESENTATION OF A POSTER AT EGS-AGU-EUG JOINT ASSEMBLY 

We have attended the EGS-AGU-EUG Joint Assembly that was held in Nice, France 
between 6 and 1 1 April 2003. A poster illustrating the structure of the code and the most 
recent result was presented. (Dalton D. Schnack and Roberto Lionello ” Developing an 
MHD Algorithm on an Unstructured Mesh of Tetrahedra”, Geophysical Research 
Abstracts, Vol. 5, 04691, 2003). The text of the abstract follows: 

MH4D (Magnetohydrodynamics on a TETRAhedral Domain) is a massively- 
parallel, device-independent numerical code for the solution of the resistive and 
viscous MHD equations on an unstructured grid of tetrahedra. The unstructured grid 
allows the resolution to be increased in the regions of physical interest. 
Consequently, MH4D can model problems with a wide range of spatial scales (e.g., 
active regions in the large scale corona). A variational formulation of the differential 
operators ensures accuracy and the preservation of the analytical properties of the 
operators ( V • B « 0) and self-adjointness of the resistive and viscous operators. The 
combined semi-implicit treatment of the waves and implicit formulation of the 
diffusive operators can accommodate the wide range of time scales present in the 
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solar corona. The capability of mesh refinement and coarsening is also included. 
MH4D is currently capable of solving the resistive diffusion equation and the 
equations of hydrodynamics. We will present preliminary results. 


2. PRESENTATION OF A POSTER AT THE SHERWOOD CONFERENCE 

We have attended the 2003 International Sherwood Theory Conference that was held 
in Corpus Christi between 28 and 30 April 2003. A poster illustrating the structure of the 
code and the most recent result was presented. (R. Lionello and D. D. Schnack " MH4D: 
an MHD Algorithm on an Unstructured Tetrahedral Grid”, 1E37). The text of the 
abstract follows: 

MH4D (Magnetohydrodynamics on a TETRAhedral Domain) is a massively- 
parallel, device-independent numerical code for the solution of the resistive and 
viscous MHD equations on an unstructured grid of tetrahedra. The unstructured grid 
allows the resolution to be increased in the regions of physical interest. 
Consequently, MH4D can model problems with a wide range of spatial scales (e.g., 
active regions in the large scale corona). A variational formulation of the differential 
operators ensures accuracy and the preservation of the analytical properties of the 
operators ( V • B « 0), and self-adjointness of the resistive and viscous operators. The 
combined semi-implicit treatment of the waves and implicit formulation of the 
diffusive operators can accommodate the wide range of time scales present in the 
solar corona. The capability of mesh refinement and coarsening is also included. 
MH4D is currently capable of solving the resistive diffusion equation and the 
equations of hydrodynamics. We are currently implementing the MHD operators and 
will present preliminary results. 


3. CODE DEVELOPMENT 

MH4D can now solve the full resistive and viscous MHD equations: 


(1) 

B * Vx A 

(2) 

J-VxB 

(3) 

dA 

— *vxB-nVxVxA 
dt 

(4) 

— + V-(pv)-0 
dt 

(5) 

f£ + y(/>v)«-(y- l)pV-v 

dt 

(6) 

dx \ 

p — + vVvl«-V/? + JxB + V-(vpVv) 


This set of equations can be advanced explicitly or semi-implicitly. We describe now 
how the magnetic operators and the semi-implicit solve have been implemented. 
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3.1 Implementation of the MHD Equations. 

We have implemented subroutines to calculate the advective term in the induction 
equation for the vector potential and the Lorentz’s force term in the momentum equation. 
The advection term, calculated on vertex i, has been written as 


( 7 ) 


(vxB), - v,x 


i». 

IUJ - 1 


where the value of the magnetic field, B/(,), is the average of all the values of the 
magnetic field located on the centroids surrounding vertex /. 

The expression for the Lorentz’s force term is 

1 Ni 

( 8 ) - 2 B m v w>- 

'/(O-l 

j) is the volume of the centroid / surrounding vertex / and V i is the volume of the dual 

mesh element centered on vertex i. This volume average has been introduced to preserve 
the self-adjointness of the numerical representation of the MHD operator. 

3.1 Implementation of the Semi-Implicit Algorithm 

In the lower corona, the Alfv&i and sound speed are large compared with the speed 
of the flow. An explicit algorithm would require a very small time-step to follow the 
propagation of short-wavelength, high frequency waves that are of minimal importance 
for the dynamics of the systems we want to model. For this reason we apply a semi- 
implicit treatment to the wave-like terms. A semi-implicit formulation avoids the time- 
step restriction imposed by high frequency compressional and torsional Alfv6n waves 
and allows for an arbitrarily large time-step with respect to the normal modes of the 
system.. In order to implement the semi-implicit method, the momentum equation has 
been modified as follows: 


(9) p— -F + aA/S'— , 

dt dt 

where F represents the explicit forces appearing in Eq. (6), a is a numerical factor, and 
S is the semi-implicit operator. Our choice for S is the following 

(10) S-v-S-VxVx(vxB)xB. + Vtf>V-v. 

S is the sum of the linearized MHD waves and sound waves operators. This choice of S 
is motivated by the fact that the Alfv6n time scale determines the stiffness of the problem 
in the lower corona. However S is easier to invert than the full MHD operator. The 
expression for S minimizes the following functional 

(11) /(v)«- [vx(vxB)| 2 + )p(V-v) 2 +2S-vJf. 

2 i/ 


We have made sure that the numerical discretization of S preserves the self-adjointess of 
the analytical operator. The operator is inverted using the same parallel, preconditioned 
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CG solver implemented to invert the resistive diffusion operator in the induction 
equation. 

4. VALIDATION OF THE NEW FEATURES IMPLEMENTED IN THE CODE 
4.1 The Sod’s Problem 


We have conducted another test of the implementation of the hydrodynamic 
equations in MH4D by reproducing the Sod’s problem (Sod, G. A. 1978, J. Comp. Phys., 
27, 1). This consists of the following: at t = 0 a diaphgram separates a tube in two regions 
filled with a fluid with different densities and pressures. 


( 12 ) 


Pleft“ 10 
Tzleft “ 0.0 


Pright * 0. 1 25 
Pright” 01 • 

. v zright=0 0 


The specific heat ratios y is set to 1 .4. At f > 0 the diaphgram is broken. We have 
designed a cylindrically shaped grid of radius 0.5 and height 1 with 10300 vertices and 
50525 cells. We have evolved the hydrodynamic equations in MH4D stopping before any 
wave has reached the left or right boundary. The results of the simulation are graphically 
presented in Fig. 1 (velocity), 2 (pressure), and 3 (density). A shock wave propagates to 
the right, a rarefaction wave to the left, and a contact discontinuity, visible in the density 
plot (Fig. 3), follows the shock. 
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FIG. 1. The component of the velocity along the tube in the Sod’s problem at 
different instants in time. After the diaphgram is broken a shockwave propagates from the 
left to the right. 
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FIG. 2. Evolution of the pressure in the Sod’s problem. A shock wave moves 
rightward and a rarefaction wave leftward. 
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FIG. 3. Evolution of the density in the Sod’s problem. We notice a shock wave 
propagating to the right followed by a contact discontinuity. A rarefaction wave moves to 
towards the left. 


4.2 Alfven Waves in a Cylinder 

We have tested the implementation of the MHD equations by reproducing 
eigenfunctions of the linearized equations (i.e. Alfven waves). The equations for linear 
Alfven waves are: 


( 13 ) 

( 14 ) 


d6 A 
dt 
ddv 

p ir 


dvxB, 

b<5JxB. 


We have used 8 to indicate the perturbed quantities. A solution (eigenvector) in a 
cylinder where B = Bi is 


( 15 ) 


dAf^ef (r) cos cot sin kz 
■ =£kf(r) cos cot cos Az, 

- e k f ( r ) sin cot sin kz 


with angular speed 
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( 16 ) 



f{r) is an arbitrary function of r. Equation (14) represents a torsional Alfv^n wave. 


For k = n, p = 1,5= 1, the theoretical period is T=2. For a cylindrical computational 
domain of radius and height equal to one, we have designed a grid of 420 vertices and 
1692 cells. We have applied the perturbation as in Eq. (15) and advanced Eqs. (1-6) with 
MH4D using both the explicit and the semi-implicit time advancement with no noticeable 
difference between the two. In Fig. 4 we show the time history of the kinetic energy, 
which shows a period of 1, as expected. In Figs. 5 and 6 we show respectively the 
poloidal component of the magnetic field at four instants in time. 



FIG. 4. Kinetic energy as a function of time for the case of the torsional Alfven wave 
eigenmode, in a cylinder. 
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FIG. 5. The poloidal magnetic field is shown at four instants in time in the 
simulation of a torsional AlfV6n wave in a cylinder obtained with MH4D. The vectors are 
colored according to their magnitude. 
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FIG. 6. The poloidal velocity field is shown at four instants in time in the simulation 
of a torsional Alfv6n wave in a cylinder obtained with MH4D. The vectors are colored 
according to their magnitude. 
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